knitr::opts_knit$set(root.dir = "/home/francescoc/Desktop/scGraphVerse/data/",message=FALSE, warning=FALSE)

Load Libraries and data

library(RColorBrewer)
library(DT)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 
## Attaching package: 'SeuratObject'
## 
## The following object is masked from 'package:DT':
## 
##     JS
## 
## The following objects are masked from 'package:base':
## 
##     intersect, t
## 
## 
## Attaching package: 'Seurat'
## 
## The following object is masked from 'package:DT':
## 
##     JS
library(STRINGdb)
library(SingleCellExperiment)
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following object is masked from 'package:SeuratObject':
## 
##     intersect
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:sp':
## 
##     %over%
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## 
## Attaching package: 'SummarizedExperiment'
## 
## The following object is masked from 'package:Seurat':
## 
##     Assays
## 
## The following object is masked from 'package:SeuratObject':
## 
##     Assays
library(scGraphVerse)
## Warning: replacing previous import 'igraph::components' by 'Seurat::components'
## when loading 'scGraphVerse'
library(BiocParallel)
library(GENIE3)

reticulate::use_python("/usr/bin/python3", required = TRUE)
arboreto <- reticulate::import("arboreto.algo")
pandas <- reticulate::import("pandas")
numpy <- reticulate::import("numpy")

time <- list()

ddir <- "/home/francescoc/Desktop/NodeVerse/analysis/GENIE3/data"
pdir <- "/home/francescoc/Desktop/NodeVerse/analysis/GENIE3/plots"

Adjacency and Count matrices

adjm <- as.matrix(read.table("./../analysis/adjm_top1200.txt"))
colnames(adjm) <- rownames(adjm)

as.data.frame(adjm) %>%
  datatable(extensions = 'Buttons',
            options = list(
              dom = 'Bfrtip',
              buttons = c('csv', 'excel'),
              scrollX = TRUE,
              pageLength = 10), 
            caption = "Ground Truth")
count_matrices <- readRDS("./../analysis/sim_n200p500.RDS")
dim(count_matrices[[1]])
## [1] 200 554
count_matrices <- lapply(count_matrices, t)

as.data.frame(count_matrices[[1]]) %>%
  slice_head(n=10) %>%
  datatable(extensions = 'Buttons',
            options = list(
              dom = 'Bfrtip',
              buttons = c('csv', 'excel'),
              scrollX = TRUE,
              pageLength = 10), 
            caption = "Simulated Count matrix")

GENIE3

Late integration

set.seed(1234)
time[["GENIE3_late_15Cores"]] <- system.time(
  genie3_late <- infer_networks(count_matrices, 
                                method="GENIE3",
                                nCores = 15)
)

genie3_late[[1]] %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 output")

Symmetrize and ROC

genie3_late_wadj <- generate_adjacency(genie3_late)
sgenie3_late_wadj <- symmetrize(genie3_late_wadj, weight_function = "mean")

as.data.frame(sgenie3_late_wadj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 symmetric weighted adjacency")
genie3_late_auc <- plotROC(sgenie3_late_wadj, adjm, plot_title = "ROC curve - GENIE3 Late Integration", is_binary = F)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

Cutoff

sgenie3_late_adj <- cutoff_adjacency(count_matrices = count_matrices,
                                     weighted_adjm_list = sgenie3_late_wadj, 
                                     n = 3,
                                     method = "GENIE3",
                                     nCores = 15)

as.data.frame(sgenie3_late_adj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 symmetric adjacency")
scores.genie3.late.all <- pscores(adjm, sgenie3_late_adj)
## Registered S3 methods overwritten by 'fmsb':
##   method    from
##   print.roc pROC
##   plot.roc  pROC

scores.genie3.late.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores late")
plotg(sgenie3_late_adj)

## TableGrob (2 x 2) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]

Consensus

consesusm <- create_consensus(sgenie3_late_adj, method="vote")
consesusu <- create_consensus(sgenie3_late_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sgenie3_late_adj, weighted_list = sgenie3_late_wadj, method = "INet", threshold = 0.05, ncores = 15)
## [1] 0.3442719
## [1] 0.06703537
scores.genie3.late <- pscores(adjm, list(consesusm,consesusu,consesunet))

scores.genie3.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote-union-iNet")

Plot comparison

compare_consensus(consensus_matrix = consesusm, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesusu, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesunet, reference_matrix = adjm, false_plot = F)

Early integration

count_matrices <- lapply(count_matrices, as.matrix)
early_matrix <- list(earlyj(count_matrices, rowg = T))

set.seed(1234)
time[["GENIE3_early_15Cores"]] <- system.time(
  genie3_early <- infer_networks(early_matrix, method="GENIE3", nCores = 15)
)

as.data.frame(genie3_early[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 early output")

Symmetrize and ROC

genie3_early_wadj <- generate_adjacency(genie3_early)
sgenie3_early_wadj <- symmetrize(genie3_early_wadj, weight_function = "mean")

as.data.frame(sgenie3_early_wadj[[1]]) %>%
  slice_head(n=30) %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 early symmetric weighted adjacency")
genie3_early_auc <- plotROC(sgenie3_early_wadj, adjm, plot_title = "ROC curve - GENIE3 Early Integration", is_binary = F)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

Cutoff

sgenie3_early_adj <- cutoff_adjacency(count_matrices = early_matrix,
                                      weighted_adjm_list = sgenie3_early_wadj, 
                                      n = 2,
                                      method = "GENIE3",
                                      nCores = 15)

as.data.frame(sgenie3_early_adj[[1]]) %>%
  slice_head(n=30) %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 early symmetric adjacency")
scores.genie3.early <- pscores(adjm, sgenie3_early_adj)

scores.genie3.early$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores early")
plotg(sgenie3_early_adj)

## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]

Plot comparison

compare_consensus(consensus_matrix = sgenie3_early_adj[[1]], reference_matrix = adjm, false_plot = F)

GRNBoost2

Late integration

set.seed(1234)
time[["grnboost_late"]] <- system.time(
  grnboost_late <- infer_networks(count_matrices, 
                                method="GRNBoost2",
                                nCores = 15)
)
## /home/francescoc/.local/lib/python3.8/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.
## Perhaps you already have a cluster running?
## Hosting the HTTP server on port 35831 instead
##   warnings.warn(
## /home/francescoc/.local/lib/python3.8/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.
## Perhaps you already have a cluster running?
## Hosting the HTTP server on port 39591 instead
##   warnings.warn(
as.data.frame(grnboost_late[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 late output")

Symmetrize and ROC

grnboost_late_wadj <- generate_adjacency(grnboost_late)
sgrnboost_late_wadj <- symmetrize(grnboost_late_wadj, weight_function = "mean")

as.data.frame(sgrnboost_late_wadj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetric weighted adjacency")
grnboost_late_auc <- plotROC(sgrnboost_late_wadj, adjm, plot_title = "ROC curve - grnboost Late Integration")
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

Cutoff

sgrnboost_late_adj <- cutoff_adjacency(count_matrices = count_matrices,
                                     weighted_adjm_list = sgrnboost_late_wadj, 
                                     n = 3,
                                     method = "GRNBoost2",
                                     nCores = 15)

as.data.frame(sgrnboost_late_adj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetric adjacency")
scores.grnboost.late.all <- pscores(adjm, sgrnboost_late_adj)

scores.grnboost.late.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores late")
plotg(sgrnboost_late_adj)

## TableGrob (2 x 2) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]

Consensus

consesusm <- create_consensus(sgrnboost_late_adj, method="vote")
consesusu <- create_consensus(sgrnboost_late_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sgrnboost_late_adj, weighted_list = sgrnboost_late_wadj, method = "INet", threshold = 0.05)
## [1] 0.2844868
## [1] 0.03984126
scores.grnboost.late <- pscores(adjm, list(consesusm,consesusu,consesunet))

scores.grnboost.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote-union-iNet")

Plot comparison

compare_consensus(consensus_matrix = consesusm, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesusu, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesunet, reference_matrix = adjm, false_plot = F)

Early integration

early_matrix <- list(earlyj(count_matrices))

set.seed(1234)
time[["grnboost_early"]] <- system.time(
  grnboost_early <- infer_networks(early_matrix, method="GRNBoost2")
)

as.data.frame(grnboost_early[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 early output")

Symmetrize and ROC

grnboost_early_wadj <- generate_adjacency(grnboost_early)
sgrnboost_early_wadj <- symmetrize(grnboost_early_wadj, weight_function = "mean")

as.data.frame(sgrnboost_early_wadj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetric weighted adjacency")
grnboost_early_auc <- plotROC(sgrnboost_early_wadj, adjm, plot_title = "ROC curve - grnboost Early Integration")
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

Cutoff

sgrnboost_early_adj <- cutoff_adjacency(count_matrices = early_matrix,
                                      weighted_adjm_list = sgrnboost_early_wadj, 
                                      n = 2,
                                      method = "GRNBoost2",
                                      nCores = 15)
as.data.frame(sgrnboost_early_adj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetric adjacency")
scores.grnboost.early <- pscores(adjm, sgrnboost_early_adj)

scores.grnboost.early$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores early")
plotg(sgrnboost_early_adj)

## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]

Plot comparison

compare_consensus(consensus_matrix = sgrnboost_early_adj[[1]], reference_matrix = adjm, false_plot = F)

Joint Integration

Joint Random Forest

#https://cran.r-project.org/src/contrib/Archive/JRF/
#install.packages("/home/francescoc/Downloads/JRF_0.1-4.tar.gz", repos = NULL, type = "source")
set.seed(1234)
time[["JRF"]] <- system.time(
  jrf_mat <- infer_networks(count_matrices, method="JRF", nCores = 15)
)
 
as.data.frame(jrf_mat[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF output")

Prepare the output

jrf_list <- list()

importance_columns <- grep("importance", names(jrf_mat[[1]]), value = TRUE)

for (i in seq_along(importance_columns)) {
# Select the 'gene1', 'gene2', and the current 'importance' column
  df <- jrf_mat[[1]][, c("gene1", "gene2", importance_columns[i])]
  
  # Rename the importance column to its original name (e.g., importance1, importance2, etc.)
  names(df)[3] <- importance_columns[i]
  
  # Add the data frame to the output list
  jrf_list[[i]] <- df
}

saveRDS(jrf_list, "./../analysis/jrf.RDS")

as.data.frame(jrf_list[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF output")

symmetrize Output and ROC

jrf_wadj <- generate_adjacency(jrf_list)
sjrf_wadj <- symmetrize(jrf_wadj, weight_function = "mean")
jrf_auc_mine <- plotROC(sjrf_wadj, adjm, plot_title = "ROC curve - JRF Late Integration", is_binary = F)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

as.data.frame(sjrf_wadj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF symmetrize output")

Generate Adjacency and Apply Cutoff

sjrf_adj <- cutoff_adjacency(count_matrices = count_matrices,
                 weighted_adjm_list = sjrf_wadj, 
                 n = 3,
                 method = "JRF")

as.data.frame(sjrf_adj[[1]]) %>%
  slice_head(n=30) %>% 
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF adjacency")

Comparison with the Ground Truth

scores.jrf.all <- pscores(adjm, sjrf_adj)

scores.jrf.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")
plotg(sjrf_adj)

## TableGrob (2 x 2) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
consesusm <- create_consensus(sjrf_adj, method="vote")
consesusu <- create_consensus(sjrf_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sjrf_adj, weighted_list = sjrf_wadj, method = "INet", threshold = 0.1, ncores = 15)
## [1] 0.1583007
## [1] 0.02881552
scores.jrf <- pscores(adjm, list(consesusm, consesusu, consesunet))

scores.jrf$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote-union-iNet")
compare_consensus(consensus_matrix = consesusm, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesusu, reference_matrix = adjm, false_plot = F)

compare_consensus(consensus_matrix = consesunet, reference_matrix = adjm, false_plot = F)

Method Comparison

time_data <- data.frame(
  Method = names(time),
  Time_in_Seconds = sapply(time, function(x) {
    if ("elapsed" %in% names(x)) x["elapsed"] else NA
  })
)

time_data$Time_in_Minutes <- as.numeric(time_data$Time_in_Seconds) / 60
time_data$Time_in_Hours <- as.numeric(time_data$Time_in_Seconds) / 3600

time_data <- time_data[order(time_data$Time_in_Hours), ]
time_data$Method <- factor(time_data$Method, levels = time_data$Method)

time_data <- time_data %>%
  mutate(Method_Group = case_when(
    grepl("GENIE3", Method) ~ "GENIE3",
    grepl("GRNBoost2", Method) ~ "GRNBoost2",
    #grepl("ZILGM", Method) ~ "ZILGM",
    grepl("JRF", Method) ~ "JRF"#,
    #grepl("PCzinb", Method) ~ "PCzinb"
  ))

method_colors <- c("GENIE3" = "darkblue",
                   "GRNBoost2" = "darkgreen",
                   #"ZILGM" = "orange",
                   "JRF" = "red"#,
                   #"PCzinb" = "purple"
                   )

time_plot <- ggplot(time_data, aes(x = Method, y = Time_in_Hours, fill = Method_Group)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = sprintf("%.1f min", Time_in_Minutes)), vjust = -0.5) +
  labs(title = "Execution Time for Each Method", y = "Time (Hours)", x = NULL) +
  scale_fill_manual(values = method_colors) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  )

print(time_plot) 

library(patchwork)
mplot <- list()

for (k in c("TPR", "FPR", "Precision", "F1", "MCC")) {
  mplot[[k]] <- data.frame(
    GENIE3_late = scores.genie3.late$Statistics[[k]][[1]],
    GENIE3_early = scores.genie3.early$Statistics[[k]],
    GRNBoost2_late = scores.grnboost.late$Statistics[[k]][[1]],
    GRNBoost2_early = scores.grnboost.early$Statistics[[k]],
    #ZILGM_late = scores.zilgm.late$Statistics[[k]],
    #ZILGM_early = scores.zilgm.early$Statistics[[k]],
    #PCzinb_late = scores.pc.late$Statistics[[k]],
    #PCzinb_early = scores.pc.early$Statistics[[k]],
    JRF = scores.jrf$Statistics[[k]][[1]]
  )
}

mplot[["AUC"]] <- data.frame(
    GENIE3_late = mean(genie3_late_auc$AUC),
    GENIE3_early = genie3_early_auc$AUC,
    GRNBoost2_late = mean(grnboost_late_auc$AUC),
    GRNBoost2_early = grnboost_early_auc$AUC,
    #ZILGM_late = zilgm_late_auc$AUC,
    #ZILGM_early = zilgm_early_auc$AUC,
    JRF = jrf_auc_mine$AUC
)

plot_data <- bind_rows(lapply(names(mplot), function(metric) {
  data.frame(
    Metric = metric,
    Method = names(mplot[[metric]]),
    Value = as.numeric(mplot[[metric]][1, ])
  )
}))

plot_data <- plot_data %>%
  mutate(Method_Group = case_when(
    grepl("GENIE3", Method) ~ "GENIE3",
    grepl("GRNBoost2", Method) ~ "GRNBoost2",
    #grepl("ZILGM", Method) ~ "ZILGM",
    grepl("JRF", Method) ~ "JRF"
  )) %>%
  mutate(Method = factor(Method, levels = c(
    "GENIE3_early", "GENIE3_late", 
    "GRNBoost2_early", "GRNBoost2_late", 
    #"ZILGM_early", "ZILGM_late",
    #"PCzinb_early", "PCzinb_late", 
    "JRF"  # Ensure JRF is last
  )))

plots <- lapply(unique(plot_data$Metric), function(metric) {
  show_x_text <- metric %in% c("Accuracy", "AUC")
  
  ggplot(plot_data %>% filter(Metric == metric), aes(x = Method, y = Value, fill = Method_Group)) +
    geom_bar(stat = "identity", position = position_dodge(width = 0.8)) +
    labs(title = metric, y = "Value", x = NULL) +
    scale_y_continuous(limits = c(0, 1)) +  
    scale_fill_manual(values = method_colors) +
    theme_minimal() +
    theme(
      axis.text.x = if (show_x_text) element_text(angle = 45, hjust = 1) else element_blank(),
      axis.title.x = element_blank(),
      legend.position = "none"  
    )
})

final_plot <- (plots[[1]] | plots[[2]]) /
              (plots[[3]] | plots[[4]]) /
              (plots[[5]] | plots[[6]]) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom")

print(final_plot)

sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=it_IT.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=it_IT.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Rome
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] patchwork_1.3.0             fmsb_0.7.6                 
##  [3] doRNG_1.8.6.1               rngtools_1.5.2             
##  [5] foreach_1.5.2               GENIE3_1.28.0              
##  [7] BiocParallel_1.40.0         scGraphVerse_0.1.0         
##  [9] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
## [11] Biobase_2.66.0              GenomicRanges_1.58.0       
## [13] GenomeInfoDb_1.42.1         IRanges_2.40.1             
## [15] S4Vectors_0.44.0            BiocGenerics_0.52.0        
## [17] MatrixGenerics_1.18.1       matrixStats_1.5.0          
## [19] STRINGdb_2.18.0             Seurat_5.2.0               
## [21] SeuratObject_5.0.2          sp_2.1-4                   
## [23] lubridate_1.9.4             forcats_1.0.0              
## [25] stringr_1.5.1               dplyr_1.1.4                
## [27] purrr_1.0.4                 readr_2.1.5                
## [29] tidyr_1.3.1                 tibble_3.2.1               
## [31] ggplot2_3.5.1               tidyverse_2.0.0            
## [33] DT_0.33                     RColorBrewer_1.1-3         
## 
## loaded via a namespace (and not attached):
##   [1] spatstat.sparse_3.1-0   bitops_1.0-9            httr_1.4.7             
##   [4] doParallel_1.0.17       numDeriv_2016.8-1.1     backports_1.5.0        
##   [7] tools_4.4.2             sctransform_0.4.1       R6_2.5.1               
##  [10] iZID_0.0.1              lazyeval_0.2.2          uwot_0.2.2             
##  [13] withr_3.0.2             gridExtra_2.3           progressr_0.15.1       
##  [16] cli_3.6.3               spatstat.explore_3.3-4  fastDummies_1.7.5      
##  [19] network_1.19.0          labeling_0.4.3          ZILGM_0.1.1            
##  [22] sass_0.4.9              spatstat.data_3.1-4     ggridges_0.5.6         
##  [25] pbapply_1.7-2           WeightSVM_1.7-16        rentrez_1.2.3          
##  [28] pscl_1.5.9              parallelly_1.41.0       plotrix_3.8-4          
##  [31] rstudioapi_0.17.1       RSQLite_2.3.9           generics_0.1.3         
##  [34] shape_1.4.6.1           crosstalk_1.2.1         gtools_3.9.5           
##  [37] ica_1.0-3               spatstat.random_3.3-2   car_3.1-3              
##  [40] Matrix_1.7-1            abind_1.4-8             lifecycle_1.0.4        
##  [43] yaml_2.3.10             carData_3.0-5           gplots_3.2.0           
##  [46] SparseArray_1.6.1       Rtsne_0.17              grid_4.4.2             
##  [49] blob_1.2.4              promises_1.3.2          crayon_1.5.3           
##  [52] miniUI_0.1.1.1          lattice_0.20-44         cowplot_1.1.3          
##  [55] pillar_1.10.1           knitr_1.49              flux_0.3-0.1           
##  [58] future.apply_1.11.3     codetools_0.2-18        glue_1.8.0             
##  [61] spatstat.univar_3.1-1   data.table_1.16.4       vctrs_0.6.5            
##  [64] png_0.1-8               spam_2.11-1             bst_0.3-24             
##  [67] gtable_0.3.6            gsubfn_0.7              cachem_1.1.0           
##  [70] xfun_0.50               S4Arrays_1.6.0          mime_0.12              
##  [73] tidygraph_1.3.1         coda_0.19-4.1           survival_3.2-11        
##  [76] datastructures_0.2.9    iterators_1.0.14        fitdistrplus_1.2-2     
##  [79] ROCR_1.0-11             learn2count_0.3.2       nlme_3.1-152           
##  [82] bit64_4.6.0-1           RcppAnnoy_0.0.22        INetTool_0.1.0         
##  [85] bslib_0.8.0             irlba_2.3.5.1           KernSmooth_2.23-20     
##  [88] rpart_4.1-15            colorspace_2.1-1        DBI_1.2.3              
##  [91] gbm_2.2.2               tidyselect_1.2.1        bit_4.5.0.1            
##  [94] mpath_0.4-2.26          compiler_4.4.2          chron_2.3-62           
##  [97] glmnet_4.1-8            graph_1.84.1            DelayedArray_0.32.0    
## [100] plotly_4.10.4           scales_1.3.0            caTools_1.18.3         
## [103] multinet_4.2.1          lmtest_0.9-40           digest_0.6.37          
## [106] goftest_1.2-3           spatstat.utils_3.1-2    rmarkdown_2.29         
## [109] XVector_0.46.0          htmltools_0.5.8.1       pkgconfig_2.0.3        
## [112] fastmap_1.2.0           rlang_1.1.5             htmlwidgets_1.6.4      
## [115] UCSC.utils_1.2.0        shiny_1.10.0            farver_2.1.2           
## [118] jquerylib_0.1.4         zoo_1.8-12              jsonlite_1.8.9         
## [121] statnet.common_4.11.0   magrittr_2.0.3          Formula_1.2-6          
## [124] GenomeInfoDbData_1.2.13 dotCall64_1.2           munsell_0.5.1          
## [127] Rcpp_1.0.14             viridis_0.6.5           proto_1.0.0            
## [130] reticulate_1.40.0       sqldf_0.4-11            pROC_1.18.5            
## [133] stringi_1.8.4           JRF_0.1-4               ggraph_2.2.1           
## [136] distributions3_0.2.2    zlibbioc_1.52.0         MASS_7.3-54            
## [139] plyr_1.8.9              parallel_4.4.2          listenv_0.9.1          
## [142] ggrepel_0.9.6           deldir_2.0-4            graphlayouts_1.2.2     
## [145] splines_4.4.2           tensor_1.5              hash_2.2.6.3           
## [148] hms_1.1.3               ggpubr_0.6.0            igraph_2.1.4           
## [151] spatstat.geom_3.3-5     ggsignif_0.6.4          RcppHNSW_0.6.0         
## [154] reshape2_1.4.4          XML_3.99-0.18           evaluate_1.0.3         
## [157] tzdb_0.4.0              tweenr_2.0.3            httpuv_1.6.15          
## [160] RANN_2.6.2              polyclip_1.10-7         future_1.34.0          
## [163] scattermore_1.2         ggforce_0.4.2           broom_1.0.7            
## [166] xtable_1.8-4            RSpectra_0.16-2         rstatix_0.7.2          
## [169] later_1.4.1             viridisLite_0.4.2       memoise_2.0.1          
## [172] cluster_2.1.2           timechange_0.3.0        globals_0.16.3